using Distributions
using StatsPlots
using StatsBase
using LaTeXStrings
using CSV
using DataFrames
using StatisticalRethinking
using StatisticalRethinking: link
using LinearAlgebra
using Logging
using Random
using Turing
# setting default attributes for plots
default(label=false)
Logging.disable_logging(Logging.Warn);
Code 4.1
n = rand(Uniform(-1, 1), 1000, 16);
pos = sum.(eachrow(n));
density(pos)
Code 4.2
prod(1 .+ rand(Uniform(0, .1), 12))
1.9303572320698072
Code 4.3
u = Uniform(0, .1)
growth = prod.(eachrow(1 .+ rand(u, 10000, 12)));
density(growth; label="density")
# overlay normal distribution
μ = mean(growth)
σ = std(growth)
plot!(Normal(μ, σ); label="Normal")
Code 4.4
big = prod.(eachrow(1 .+ rand(Uniform(0, 0.5), 10000, 12)));
small = prod.(eachrow(1 .+ rand(Uniform(0, 0.01), 10000, 12)));
density([big, small]; layout=(2, 1))
Code 4.5
density(log.(big))
Code 4.6
w = 6
n = 9
p_grid = range(0, 1; length=100)
bin_dens = [pdf(Binomial(n, p), w) for p in p_grid]
uni_dens = [pdf(Uniform(0, 1), p) for p in p_grid];
posterior = bin_dens .* uni_dens
posterior /= sum(posterior);
Code 4.7
d = DataFrame(CSV.File("data/Howell1.csv"));
Code 4.8
describe(d)
4 rows × 7 columns
| variable | mean | min | median | max | nmissing | eltype | |
|---|---|---|---|---|---|---|---|
| Symbol | Float64 | Real | Float64 | Real | Int64 | DataType | |
| 1 | height | 138.264 | 53.975 | 148.59 | 179.07 | 0 | Float64 |
| 2 | weight | 35.6106 | 4.25242 | 40.0578 | 62.9926 | 0 | Float64 |
| 3 | age | 29.3444 | 0.0 | 27.0 | 88.0 | 0 | Float64 |
| 4 | male | 0.472426 | 0 | 0.0 | 1 | 0 | Int64 |
Code 4.9
precis(d)
┌────────┬────────────────────────────────────────────────────────────┐ │ param │ mean std 5.5% 50% 94.5% histogram │ ├────────┼────────────────────────────────────────────────────────────┤ │ height │ 138.264 27.6024 81.1086 148.59 165.735 ▁▁▁▂▂▂▂▂▂██▆▁ │ │ weight │ 35.6106 14.7192 9.3607 40.0578 54.5029 ▁▃▄▄▃▂▃▆██▅▃▁ │ │ age │ 29.3444 20.7469 1.0 27.0 66.135 █▆▆▆▆▃▃▁▁ │ │ male │ 0.4724 0.4997 0.0 0.0 1.0 █▁▁▁▁▁▁▁▁▁█ │ └────────┴────────────────────────────────────────────────────────────┘
Code 4.10
d.height
544-element Vector{Float64}:
151.765
139.7
136.525
156.845
145.415
163.83
149.225
168.91
147.955
165.1
154.305
151.13
144.78
⋮
156.21
152.4
162.56
114.935
67.945
142.875
76.835
145.415
162.56
156.21
71.12
158.75
Code 4.11
d2 = d[d.age .>= 18,:];
Code 4.12
plot(Normal(178, 20); xlim=(100, 250))
Code 4.13
plot(Uniform(0, 50), xlim=(-10, 60), ylim=(0, 0.03))
Code 4.14
size = 10_000
sample_μ = rand(Normal(178, 20), size)
sample_σ = rand(Uniform(0, 50), size);
prior_h = [rand(Normal(μ, σ)) for (μ, σ) in zip(sample_μ, sample_σ)];
p1 = density(sample_μ; title="μ")
p2 = density(sample_σ; title="σ")
p3 = density(prior_h; title="prior_h")
plot(p1, p2, p3, layout=@layout [a b; c])
Code 4.15
sample_μ = rand(Normal(178, 100), size)
prior_h = [rand(Normal(μ, σ)) for (μ, σ) in zip(sample_μ, sample_σ)];
density(prior_h)
vline!([0, 272])
Code 4.16
μ_list = range(150, 160; length=100)
σ_list = range(7, 9; length=100)
log_likelihood = [
sum(logpdf(Normal(μ, σ), d2.height))
for μ ∈ μ_list, σ ∈ σ_list
]
log_prod = log_likelihood .+ [
logpdf(Normal(178, 20), μ) + logpdf(Uniform(0, 50), σ)
for μ ∈ μ_list, σ ∈ σ_list
];
max_prod = maximum(log_prod)
prob = @. exp(log_prod - max_prod);
Code 4.17
# note the transposition, that's due to Julia matrix order
contour(μ_list, σ_list, prob')
Code 4.18
heatmap(μ_list, σ_list, prob')
Code 4.19
indices = collect(Iterators.product(1:length(μ_list), 1:length(σ_list)));
sample_idx = wsample(vec(indices), vec(prob), 10_000; replace=true)
sample_μ = μ_list[first.(sample_idx)]
sample_σ = σ_list[last.(sample_idx)];
Code 4.20
scatter(sample_μ, sample_σ; alpha=0.1)
Code 4.21
p1 = density(sample_μ)
p2 = density(sample_σ)
plot(p1, p2, layout=(2,1))
Code 4.22
println(hpdi(sample_μ, alpha=0.11))
println(hpdi(sample_σ, alpha=0.11))
[153.83838383838383, 155.15151515151516] [7.262626262626263, 8.191919191919192]
Code 4.23
d3 = sample(d2.height, 20);
Code 4.24
μ_list = range(150, 170; length=200)
σ_list = range(4, 20; length=200)
log_likelihood = [
sum(logpdf(Normal(μ, σ), d3))
for μ ∈ μ_list, σ ∈ σ_list
]
log_prod = log_likelihood .+ [
logpdf(Normal(178, 20), μ) + logpdf(Uniform(0, 50), σ)
for μ ∈ μ_list, σ ∈ σ_list
]
max_prod = maximum(log_prod)
prob2 = @. exp(log_prod - max_prod)
indices = collect(Iterators.product(1:length(μ_list), 1:length(σ_list)));
sample2_idx = wsample(vec(indices), vec(prob2), 10_000; replace=true)
sample2_μ = μ_list[first.(sample2_idx)]
sample2_σ = σ_list[last.(sample2_idx)]
scatter(sample2_μ, sample2_σ; alpha=0.1)
Code 4.25
density(sample2_σ)
μ = mean(sample2_σ)
σ = std(sample2_σ)
plot!(Normal(μ, σ); label="Normal")
Code 4.26
d = DataFrame(CSV.File("data/Howell1.csv"));
d2 = d[d.age .>= 18,:];
Code 4.27
@model function model_height(height)
μ ~ Normal(178, 20)
σ ~ Uniform(0, 50)
height ~ Normal(μ, σ)
end
model_height (generic function with 2 methods)
Code 4.28
m4_1 = sample(model_height(d2.height), NUTS(), 1000)
Chains MCMC chain (1000×14×1 Array{Float64, 3}):
Iterations = 501:1:1500
Number of chains = 1
Samples per chain = 1000
Wall duration = 7.07 seconds
Compute duration = 7.07 seconds
parameters = μ, σ
internals = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size
Summary Statistics
parameters mean std naive_se mcse ess rhat ⋯
Symbol Float64 Float64 Float64 Float64 Float64 Float64 ⋯
μ 154.6092 0.3988 0.0126 0.0113 1239.9501 0.9990 ⋯
σ 7.7657 0.3138 0.0099 0.0092 1149.2957 0.9995 ⋯
1 column omitted
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
μ 153.8444 154.3255 154.6183 154.8714 155.4424
σ 7.1909 7.5551 7.7388 7.9676 8.4361
Code 4.29
display.(describe(m4_1; q=[0.055, 0.945]));
Summary Statistics parameters mean std naive_se mcse ess rhat ⋯ Symbol Float64 Float64 Float64 Float64 Float64 Float64 ⋯ μ 154.6092 0.3988 0.0126 0.0113 1239.9501 0.9990 ⋯ σ 7.7657 0.3138 0.0099 0.0092 1149.2957 0.9995 ⋯ 1 column omitted
Quantiles parameters 5.5% 94.5% Symbol Float64 Float64 μ 153.9801 155.2242 σ 7.2873 8.2936
Code 4.30
init_vals = [mean(d2.height), std(d2.height)]
chain = sample(model_height(d2.height), NUTS(), 1000, init_theta=init_vals)
Chains MCMC chain (1000×14×1 Array{Float64, 3}):
Iterations = 501:1:1500
Number of chains = 1
Samples per chain = 1000
Wall duration = 0.65 seconds
Compute duration = 0.65 seconds
parameters = μ, σ
internals = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size
Summary Statistics
parameters mean std naive_se mcse ess rhat ⋯
Symbol Float64 Float64 Float64 Float64 Float64 Float64 ⋯
μ 154.6284 0.4095 0.0129 0.0129 813.0170 0.9999 ⋯
σ 7.7597 0.2926 0.0093 0.0089 984.3953 0.9992 ⋯
1 column omitted
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
μ 153.8513 154.3544 154.6117 154.9108 155.4373
σ 7.2344 7.5607 7.7465 7.9384 8.3772
Code 4.31
@model function model_height(height)
μ ~ Normal(178, 0.1)
σ ~ Uniform(0, 50)
height ~ Normal(μ, σ)
end
m4_2 = sample(model_height(d2.height), NUTS(), 1000)
display.(describe(m4_2; q=[0.055, 0.945]));
Summary Statistics parameters mean std naive_se mcse ess rhat ⋯ Symbol Float64 Float64 Float64 Float64 Float64 Float64 ⋯ μ 177.8637 0.1023 0.0032 0.0030 811.9349 1.0033 ⋯ σ 24.5429 0.9562 0.0302 0.0276 834.7287 0.9991 ⋯ 1 column omitted
Quantiles parameters 5.5% 94.5% Symbol Float64 Float64 μ 177.6973 178.0280 σ 23.0719 26.0606
Code 4.32
cov(hcat(m4_1[:μ], m4_1[:σ]))
2×2 Matrix{Float64}:
0.159072 0.00303165
0.00303165 0.09848
Code 4.33
c = cov(hcat(m4_1[:μ], m4_1[:σ]))
cov2cor(c, diag(c))
2×2 Matrix{Float64}:
1.0 0.193526
0.193526 1.0
Code 4.34
# resetrange is needed due to bug in MCMCChains: https://github.com/TuringLang/MCMCChains.jl/issues/324
# once it will be fixed, direct sampling from the chain will be possible
samp_chain = sample(resetrange(m4_1), 10_000)
samp_df = DataFrame(samp_chain)
first(samp_df, 5)
5 rows × 2 columns
| μ | σ | |
|---|---|---|
| Float64 | Float64 | |
| 1 | 154.911 | 7.40381 |
| 2 | 155.033 | 7.69794 |
| 3 | 154.885 | 8.00863 |
| 4 | 154.094 | 7.39209 |
| 5 | 154.939 | 7.69093 |
Code 4.35
precis(samp_df)
┌───────┬─────────────────────────────────────────────────────────┐ │ param │ mean std 5.5% 50% 94.5% histogram │ ├───────┼─────────────────────────────────────────────────────────┤ │ μ │ 154.611 0.4013 153.968 154.62 155.247 ▁▁▂▄▇▇█▇▄▂▁▁ │ │ σ │ 7.7652 0.3129 7.2866 7.7433 8.2927 ▁▁▃▆█▆▅▂▁▁ │ └───────┴─────────────────────────────────────────────────────────┘
Code 4.36
data = hcat(m4_1[:μ], m4_1[:σ])
μ = mean(data, dims=1)
σ = cov(data)
mvn = MvNormal(vec(μ), σ)
post = rand(mvn, 10_000);
print(mvn)
FullNormal( dim: 2 μ: [154.6091697298638, 7.765669889245865] Σ: [0.15907154466725013 0.003031652659777947; 0.003031652659777947 0.09848001526298829] )
Code 4.37
d = DataFrame(CSV.File("data/Howell1.csv"));
d2 = d[d.age .>= 18,:];
# fancy way of doing scatter(d2.weight, d2.height)
@df d2 scatter(:weight, :height)
Code 4.38
Random.seed!(2971)
N = 100
a = rand(Normal(178, 20), N)
b = rand(Normal(0, 10), N);
Code 4.39
p = hline([0, 272]; ylims=(-100, 400), xlabel="weight", ylabel="hegiht")
title!(L"\beta \sim \mathcal{N}(\mu=0,\sigma=10)")
x_mean = mean(d2.weight)
xlims = extrema(d2.weight) # getting min and max in one pass
for (α, β) ∈ zip(a, b)
plot!(x -> α + β * (x - x_mean); xlims=xlims, c=:black, alpha=0.3)
end
display(p)
Code 4.40
b = rand(LogNormal(0, 1), 10_000)
density(b, xlims=(0, 5), bandwidth=0.1)
Code 4.41
Random.seed!(2971)
N = 100
a = rand(Normal(178, 20), N)
b = rand(LogNormal(0, 1), N);
Code 4.42
d = DataFrame(CSV.File("data/Howell1.csv"));
d2 = d[d.age .>= 18,:]
xbar = mean(d2.weight)
@model function height_regr_model(weight, height)
a ~ Normal(178, 20)
b ~ LogNormal(0, 1)
μ = @. a + b * (weight - xbar)
σ ~ Uniform(0, 50)
height ~ MvNormal(μ, σ)
end
m4_3 = sample(height_regr_model(d2.weight, d2.height), NUTS(), 1000)
m4_3 = resetrange(m4_3);
Code 4.43
@model function height_regr_model_exp(weight, height)
a ~ Normal(178, 20)
log_b ~ Normal(0, 1)
μ = @. a + exp(log_b) * (weight - xbar)
σ ~ Uniform(0, 50)
height ~ MvNormal(μ, σ)
end
m4_3b = sample(height_regr_model_exp(d2.weight, d2.height), NUTS(), 1000);
Code 4.44
m4_3_df = DataFrame(m4_3)
precis(m4_3_df)
┌───────┬────────────────────────────────────────────────────────┐ │ param │ mean std 5.5% 50% 94.5% histogram │ ├───────┼────────────────────────────────────────────────────────┤ │ a │ 154.606 0.2625 154.175 154.616 155.014 ▁▁▂▅██▅▂▁▁ │ │ b │ 0.9014 0.0423 0.835 0.9021 0.9665 ▁▃██▃▁ │ │ σ │ 5.1071 0.1978 4.8019 5.0921 5.4553 ▁▁▅█▅▂▁ │ └───────┴────────────────────────────────────────────────────────┘
Code 4.45
round.(cov(Matrix(m4_3_df)), digits=3)
3×3 Matrix{Float64}:
0.069 -0.0 -0.001
-0.0 0.002 -0.0
-0.001 -0.0 0.039
Code 4.46
p = @df d2 scatter(:weight, :height; alpha=0.3)
chain = resetrange(m4_3)
samples = sample(chain, 1000)
a_map = mean(samples[:a])
b_map = mean(samples[:b])
plot!(x -> a_map + b_map*(x-xbar))
Code 4.47
post = sample(m4_3, 1000)
post_df = DataFrame(post)
post_df[1:5,:]
5 rows × 3 columns
| a | b | σ | |
|---|---|---|---|
| Float64 | Float64 | Float64 | |
| 1 | 155.241 | 0.939497 | 5.32578 |
| 2 | 154.435 | 0.890103 | 5.15068 |
| 3 | 154.393 | 0.95766 | 5.12624 |
| 4 | 154.644 | 0.975131 | 5.02818 |
| 5 | 154.719 | 0.833828 | 5.12795 |
Code 4.48
N = 10
dN = d2[1:N,:]
@model function height_regr_model_N(weight, height)
a ~ Normal(178, 20)
b ~ LogNormal(0, 1)
m_weight = mean(weight)
μ = @. a + b * (weight - m_weight)
σ ~ Uniform(0, 50)
height ~ MvNormal(μ, σ)
end
mN = sample(height_regr_model(dN.weight, dN.height), NUTS(), 1000)
mN = resetrange(mN);
Code 4.49
post = sample(mN, 20)
post_df = DataFrame(post);
xlims = extrema(d2.weight)
ylims = extrema(d2.height)
p = @df dN scatter(:weight, :height; xlims=xlims, ylims=ylims)
title!("N = $N"; xlab="weight", ylab="height")
x_mean = mean(dN.weight)
for (a, b) ∈ zip(post_df.a, post_df.b)
plot!(x -> a + b * (x-x_mean); c="black", alpha=0.3)
end
display(p)
Code 4.50
post = sample(m4_3, 1000)
post_df = DataFrame(post)
μ_at_50 = @. post_df.a + post_df.b * (50 - xbar);
Code 4.51
density(μ_at_50; lw=2, xlab=L"\mu|weight=50")
Code 4.52
PI(μ_at_50)
2-element Vector{Float64}:
158.5792439935818
159.65866731906047
Code 4.53
μ = link(post_df, [:a :b], d2.weight, xbar);
μ = hcat(μ...);
Base.size(μ), μ[1:5,1]
((1000, 352), [157.56293534189655, 157.70353697789378, 156.97988672219105, 156.85965609569527, 157.2632797224988])
Code 4.54
weight_seq = 25:70
μ = link(post_df, [:a :b], weight_seq, xbar);
μ = hcat(μ...);
Base.size(μ), μ[1:5,1]
((1000, 46), [137.04072934028633, 136.64103095564798, 135.80904017330352, 136.73820023564727, 137.81265964505803])
Code 4.55
p = plot()
for i in 1:100
scatter!(weight_seq, μ[i,:]; c=:blue, alpha=0.2)
end
display(p)
Code 4.56
μ_mean = mean.(eachcol(μ))
μ_PI = PI.(eachcol(μ))
μ_PI = vcat(μ_PI'...);
Code 4.57
@df d2 scatter(:weight, :height; alpha=0.2, xlab="weight", ylab="height")
plot!(weight_seq, [μ_mean μ_mean]; c=:black, fillrange=μ_PI, fillalpha=0.3)
Code 4.58
post = sample(m4_3, 1000)
post = DataFrame(post)
weight_seq = 25:70
μ = map(w -> post.a + post.b * (w - xbar), weight_seq)
μ = hcat(μ...)
μ_mean = mean.(eachcol(μ))
μ_CI = PI.(eachcol(μ));
Code 4.59
sim_height = simulate(post, [:a, :b, :σ], weight_seq .- xbar);
Base.size(sim_height), sim_height[1:5,1]
((1000, 46), [131.0600330744535, 131.07613382456364, 144.0054107126664, 133.0809884065339, 139.1993923582497])
Code 4.60
height_PI = PI.(eachcol(sim_height))
height_PI = vcat(height_PI'...);
Code 4.61
@df d2 scatter(:weight, :height; alpha=0.2, xlab="weight", ylab="height")
plot!(weight_seq, [μ_mean μ_mean]; c=:black, fillrange=μ_PI, fillalpha=0.3)
plot!(weight_seq, [μ_mean μ_mean]; c=:black, fillrange=height_PI, fillalpha=0.3)
Code 4.62
post = sample(m4_3, 10_000)
post = DataFrame(post)
sim_height = simulate(post, [:a, :b, :σ], weight_seq .- xbar)
height_PI = PI.(eachcol(sim_height))
height_PI = vcat(height_PI'...);
@df d2 scatter(:weight, :height; alpha=0.2, xlab="weight", ylab="height")
plot!(weight_seq, [μ_mean μ_mean]; c=:black, fillrange=μ_PI, fillalpha=0.3)
plot!(weight_seq, [μ_mean μ_mean]; c=:black, fillrange=height_PI, fillalpha=0.3)
Code 4.63
post = sample(m4_3, 1000)
post = DataFrame(post)
sim_height = [
[
rand(Normal(a + b * (w - xbar), σ))
for (a, b, σ) ∈ zip(post.a, post.b, post.σ)
]
for w ∈ weight_seq
]
sim_height = hcat(sim_height...)
height_PI = PI.(eachcol(sim_height));
height_PI = vcat(height_PI'...);
@df d2 scatter(:weight, :height; alpha=0.2, xlab="weight", ylab="height")
plot!(weight_seq, [μ_mean μ_mean]; c=:black, fillrange=μ_PI, fillalpha=0.3)
plot!(weight_seq, [μ_mean μ_mean]; c=:black, fillrange=height_PI, fillalpha=0.3)
Code 4.64
d = DataFrame(CSV.File("data/Howell1.csv"))
scatter(d.weight, d.height; alpha=0.3)
Code 4.65
d[!, :weight_s] = standardize(ZScoreTransform, d.weight)
d[!, :weight_s2] = d.weight_s.^2;
@model function height_regr_m2(weight_s, weight_s2, height)
a ~ Normal(178, 20)
b1 ~ LogNormal(0, 1)
b2 ~ Normal(0, 1)
μ = @. a + b1 * weight_s + b2 * weight_s2
σ ~ Uniform(0, 50)
height ~ MvNormal(μ, σ)
end
m4_5 = sample(height_regr_m2(d.weight_s, d.weight_s2, d.height), NUTS(), 1000)
m4_5 = resetrange(m4_5)
m4_5_df = DataFrame(m4_5);
Code 4.66
precis(m4_5_df)
┌───────┬───────────────────────────────────────────────────────────┐ │ param │ mean std 5.5% 50% 94.5% histogram │ ├───────┼───────────────────────────────────────────────────────────┤ │ a │ 146.054 0.3497 145.501 146.059 146.576 ▁▁▂▃▅██▇▄▂▁▁▁ │ │ b1 │ 21.7228 0.2941 21.2498 21.7259 22.1891 ▁▁▁▄▇█▇▄▂▁▁ │ │ b2 │ -7.8033 0.2589 -8.1943 -7.8236 -7.3754 ▁▁▂▄█▆▅▂▁ │ │ σ │ 5.8137 0.179 5.5276 5.8055 6.1068 ▁▁▁▃▆██▆▄▂▁▁ │ └───────┴───────────────────────────────────────────────────────────┘
Code 4.67
Random.seed!(1)
df = sample(m4_5_df, 1000)
weight_seq = range(-2.2, 2; length=30)
mu = link(df, (r, x) -> r.a + r.b1*x + r.b2*x^2, weight_seq)
mu = hcat(mu...)
mu_mean = mean.(eachcol(mu))
mu_PI = PI.(eachcol(mu))
mu_PI = vcat(mu_PI'...)
sim_height = simulate(df, (r, x) -> Normal(r.a + r.b1*x + r.b2*x^2, r.σ), weight_seq)
sim_height = vcat(sim_height'...);
height_PI = PI.(eachcol(sim_height))
height_PI = vcat(height_PI'...);
Code 4.68
p_square = @df d scatter(:weight_s, :height; alpha=0.3, title="Square poly")
plot!(weight_seq, mu_mean; c=:black)
plot!(weight_seq, [mu_mean mu_mean]; c=:black, fillrange=mu_PI, fillalpha=0.3)
plot!(weight_seq, [mu_mean mu_mean]; c=:black, fillrange=height_PI, fillalpha=0.3)
Code 4.69
d[!, :weight_s3] = d.weight_s.^3;
@model function height_regr_m3(weight_s, weight_s2, weight_s3, height)
a ~ Normal(178, 20)
b1 ~ LogNormal(0, 1)
b2 ~ Normal(0, 1)
b3 ~ Normal(0, 1)
μ = @. a + b1 * weight_s + b2 * weight_s2 + b3 * weight_s3
σ ~ Uniform(0, 50)
height ~ MvNormal(μ, σ)
end
m4_6 = sample(height_regr_m3(d.weight_s, d.weight_s2, d.weight_s3, d.height), NUTS(), 1000)
m4_6 = resetrange(m4_6)
m4_6_df = DataFrame(m4_6)
precis(m4_6_df)
┌───────┬───────────────────────────────────────────────────────────┐ │ param │ mean std 5.5% 50% 94.5% histogram │ ├───────┼───────────────────────────────────────────────────────────┤ │ a │ 146.399 0.3158 145.888 146.405 146.898 ▁▁▁▃▅██▆▃▁▁ │ │ b1 │ 15.2346 0.4765 14.4594 15.2243 15.9857 ▁▁▁▄█▅▁▁ │ │ b2 │ -6.213 0.2561 -6.6352 -6.2127 -5.8004 ▁▁▂▅██▄▂▁▁ │ │ b3 │ 3.5764 0.2239 3.2084 3.584 3.9364 ▁▂▄██▃▁▁ │ │ σ │ 4.8605 0.147 4.6456 4.8546 5.0868 ▁▁▁▄▇█▇▄▂▁▁▁▁ │ └───────┴───────────────────────────────────────────────────────────┘
Random.seed!(1)
df = sample(m4_6_df, 1000)
weight_seq = range(-2.2, 2; length=30)
func = (r, x) -> r.a + r.b1*x + r.b2*x^2 + r.b3*x^3
mu = link(df, func, weight_seq)
mu = hcat(mu...)
mu_mean = mean.(eachcol(mu))
mu_PI = PI.(eachcol(mu))
mu_PI = vcat(mu_PI'...)
func = (r, x) -> Normal(r.a + r.b1*x + r.b2*x^2 + r.b3*x^3, r.σ)
sim_height = simulate(df, func, weight_seq)
sim_height = vcat(sim_height'...);
height_PI = PI.(eachcol(sim_height))
height_PI = vcat(height_PI'...);
p_cubic = @df d scatter(:weight_s, :height; alpha=0.3, title="Cubic poly")
plot!(weight_seq, mu_mean; c=:black)
plot!(weight_seq, [mu_mean mu_mean]; c=:black, fillrange=mu_PI, fillalpha=0.3)
plot!(weight_seq, [mu_mean mu_mean]; c=:black, fillrange=height_PI, fillalpha=0.3)
plot(p_square, p_cubic; layout=(1, 2))
Code 4.70 and 4.71
Looks like Julia plots don't support change of ticks proposed in the book. But much more natural way will be to remap values we're plotting back to the original scale. Example of this is below.
μ = mean(d.weight)
σ = std(d.weight)
w = @. d.weight_s * σ + μ
scatter(w, d.height; alpha=0.3)
w_s = @. weight_seq * σ + μ
plot!(w_s, mu_mean; c=:black)
plot!(w_s, [mu_mean mu_mean]; c=:black, fillrange=mu_PI, fillalpha=0.3)
plot!(w_s, [mu_mean mu_mean]; c=:black, fillrange=height_PI, fillalpha=0.3)
Code 4.72
d = DataFrame(CSV.File("data/cherry_blossoms.csv", missingstring="NA"))
precis(d)
┌────────────┬─────────────────────────────────────────────────────────┐ │ param │ mean std 5.5% 50% 94.5% histogram │ ├────────────┼─────────────────────────────────────────────────────────┤ │ year │ 1408.0 350.885 867.77 1408.0 1948.23 ████████████▂ │ │ doy │ 104.54 6.407 94.43 105.0 115.0 ▁▂▆██▅▂▁ │ │ temp │ 6.1419 0.6636 5.15 6.1 7.2947 ▁▄▆█▄▂▁▁ │ │ temp_upper │ 7.1852 0.9929 5.8977 7.04 8.9023 ▂██▃▁▁▁▁ │ │ temp_lower │ 5.0989 0.8503 3.7876 5.145 6.37 ▁▁▁▂▇█▂▁ │ └────────────┴─────────────────────────────────────────────────────────┘
Code 4.73
d2 = d[completecases(d[!,[:doy]]),:]
d2 = disallowmissing(d2[!,[:year,:doy]])
num_knots = 15
knots_list = quantile(d2.year, range(0, 1; length=num_knots));
Code 4.74
using BSplines
basis = BSplineBasis(3, knots_list)
16-element BSplineBasis{Vector{Float64}}:
order: 3
breakpoints: [812.0, 1036.0, 1174.0, 1269.0, 1377.0, 1454.0, 1518.0, 1583.0, 1650.0, 1714.0, 1774.0, 1833.0, 1893.0, 1956.0, 2015.0]
Code 4.75
p1 = plot(basis)
scatter!(knots_list, repeat([1], num_knots); xlab="year", ylab="basis", legend=false)
Code 4.76
This way of calucalting bsplines is slightly slower, than shown in the book (with pre-calculated matrix) but it is much cleaner in my perspective.
You can do comparison yourself by precalculating spline matrix outside of the model and do matrix multiplication in the model instead of spline evialutaion. Example of doing this is at code block 4.79
@model function model_splines(year, doy)
w ~ MvNormal(zeros(length(basis)), 1)
a ~ Normal(100, 10)
s = Spline(basis, w)
μ = a .+ s.(year)
σ ~ Exponential(1)
doy ~ MvNormal(μ, σ)
end
m4_7 = sample(model_splines(d2.year, d2.doy), NUTS(0.65; init_ϵ = 9.765625e-5), 1000)
Chains MCMC chain (1000×30×1 Array{Float64, 3}):
Iterations = 501:1:1500
Number of chains = 1
Samples per chain = 1000
Wall duration = 12.52 seconds
Compute duration = 12.52 seconds
parameters = w[1], w[2], w[3], w[4], w[5], w[6], w[7], w[8], w[9], w[10], w[11], w[12], w[13], w[14], w[15], w[16], a, σ
internals = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size
Summary Statistics
parameters mean std naive_se mcse ess rhat ⋯
Symbol Float64 Float64 Float64 Float64 Float64 Float64 ⋯
w[1] -0.6812 0.9172 0.0290 0.0163 1813.8789 0.9991 ⋯
w[2] -0.8555 0.8626 0.0273 0.0187 2705.8346 0.9990 ⋯
w[3] 0.8252 0.8027 0.0254 0.0213 1626.8872 1.0002 ⋯
w[4] 0.6955 0.7619 0.0241 0.0167 1808.3373 0.9990 ⋯
w[5] 0.5630 0.7818 0.0247 0.0127 2370.8672 1.0008 ⋯
w[6] -1.2257 0.7915 0.0250 0.0198 2674.2308 0.9994 ⋯
w[7] 0.2788 0.7483 0.0237 0.0189 2438.4965 0.9996 ⋯
w[8] 1.8091 0.7712 0.0244 0.0163 2157.9901 1.0000 ⋯
w[9] -0.0526 0.7338 0.0232 0.0121 2400.1936 0.9991 ⋯
w[10] 1.9439 0.8019 0.0254 0.0152 2163.3095 0.9995 ⋯
w[11] 0.6797 0.7354 0.0233 0.0197 1598.7700 1.0000 ⋯
w[12] 1.2764 0.7565 0.0239 0.0146 2692.3364 0.9990 ⋯
w[13] 1.2634 0.7793 0.0246 0.0190 2460.8566 0.9994 ⋯
w[14] -0.7788 0.8162 0.0258 0.0157 2129.6796 0.9998 ⋯
w[15] -2.9017 0.8942 0.0283 0.0180 2867.7978 0.9991 ⋯
w[16] -2.7443 0.8905 0.0282 0.0188 2674.3391 0.9994 ⋯
a 104.2685 0.3292 0.0104 0.0104 1152.1117 1.0032 ⋯
σ 6.1107 0.1658 0.0052 0.0035 2781.4362 0.9991 ⋯
1 column omitted
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
w[1] -2.5031 -1.2788 -0.6836 -0.0650 1.0794
w[2] -2.5142 -1.4874 -0.8544 -0.2671 0.8429
w[3] -0.7560 0.2729 0.8494 1.3813 2.3222
w[4] -0.7929 0.1645 0.6935 1.1964 2.1724
w[5] -0.9255 0.0087 0.5661 1.0865 2.1346
w[6] -2.8200 -1.7612 -1.2422 -0.6835 0.3608
w[7] -1.2176 -0.2135 0.2979 0.7884 1.7322
w[8] 0.2887 1.2753 1.8034 2.3243 3.3166
w[9] -1.4817 -0.5690 -0.0725 0.4282 1.3626
w[10] 0.3610 1.4107 1.9837 2.5043 3.4097
w[11] -0.8514 0.1983 0.6660 1.1726 2.1750
w[12] -0.1970 0.7800 1.2584 1.8213 2.7546
w[13] -0.3652 0.7719 1.2357 1.8412 2.7361
w[14] -2.3775 -1.3282 -0.7610 -0.2397 0.8342
w[15] -4.7354 -3.4840 -2.8894 -2.3133 -1.1709
w[16] -4.5349 -3.3456 -2.7131 -2.1420 -1.0208
a 103.6152 104.0409 104.2619 104.4908 104.8871
σ 5.7975 5.9989 6.1047 6.2182 6.4409
Code 4.77
post = DataFrame(m4_7)
# convert columns w[*] into single column w
w_df = DataFrames.select(post, r"w")
post = DataFrames.select(post, Not(r"w"))
post[!,:w] = Vector.(eachrow(w_df))
# vector of 16 average w values
w_mean = mean.(eachcol(w_df))
p2 = plot(basis .* w_mean)
scatter!(knots_list, repeat([1], num_knots); xlab="year", ylab="basis × weight")
Code 4.78
μ = link(post, (r, x) -> r.a + Spline(basis, r.w)(x), d2.year)
μ = vcat(μ'...);
μ_PI = PI.(eachrow(μ))
μ_PI = vcat(μ_PI'...);
p3 = @df d2 scatter(:year, :doy; alpha=0.3)
μ_mean = mean.(eachrow(μ_PI))
plot!(d2.year, [μ_mean, μ_mean]; c=:black, fillrange=μ_PI, fillalpha=0.3, alpha=0)
plot(p1, p2, p3; layout=(3, 1))
Code 4.79
How to build the model with explicit spline matrix calculation
basis = BSplineBasis(3, knots_list)
# list of splines with 1 only at corresponding basis index
splines = [
Spline(basis, [float(idx == knot) for idx ∈ 1:length(basis)])
for knot ∈ 1:length(basis)
]
# calculate each spline for every year. Resulting matrix B is 827x16
B = [
map(s -> s(year), splines)
for year in d2.year
]
B = vcat(B'...);
# do not need years parameter anymore, all the information is in B matrix
@model function model_splines_matrix(doy)
w ~ MvNormal(zeros(length(basis)), 1)
a ~ Normal(100, 10)
μ = a .+ B * w
σ ~ Exponential(1)
doy ~ MvNormal(μ, σ)
end
m4_7alt = sample(model_splines_matrix(d2.doy), NUTS(0.65; init_ϵ = 0.0001953125), 1000)
Chains MCMC chain (1000×30×1 Array{Float64, 3}):
Iterations = 501:1:1500
Number of chains = 1
Samples per chain = 1000
Wall duration = 4.89 seconds
Compute duration = 4.89 seconds
parameters = w[1], w[2], w[3], w[4], w[5], w[6], w[7], w[8], w[9], w[10], w[11], w[12], w[13], w[14], w[15], w[16], a, σ
internals = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size
Summary Statistics
parameters mean std naive_se mcse ess rhat ⋯
Symbol Float64 Float64 Float64 Float64 Float64 Float64 ⋯
w[1] -0.6977 0.9024 0.0285 0.0185 2200.5259 0.9990 ⋯
w[2] -0.8587 0.8423 0.0266 0.0170 2112.9445 1.0009 ⋯
w[3] 0.8351 0.7801 0.0247 0.0161 2085.0565 0.9990 ⋯
w[4] 0.6791 0.6900 0.0218 0.0167 1563.5078 0.9990 ⋯
w[5] 0.5912 0.8079 0.0255 0.0170 2520.6431 0.9990 ⋯
w[6] -1.2373 0.7507 0.0237 0.0151 2291.9338 0.9995 ⋯
w[7] 0.2550 0.8036 0.0254 0.0215 1958.5399 0.9991 ⋯
w[8] 1.7987 0.8424 0.0266 0.0168 1915.4226 0.9994 ⋯
w[9] -0.0478 0.7916 0.0250 0.0161 2057.6024 0.9993 ⋯
w[10] 1.9124 0.7558 0.0239 0.0170 1717.8319 1.0001 ⋯
w[11] 0.6770 0.8081 0.0256 0.0220 1902.8984 0.9991 ⋯
w[12] 1.2699 0.7636 0.0241 0.0184 1858.3391 0.9991 ⋯
w[13] 1.2766 0.7838 0.0248 0.0145 2455.5405 0.9990 ⋯
w[14] -0.7512 0.7893 0.0250 0.0141 1702.2662 0.9990 ⋯
w[15] -2.8963 0.8600 0.0272 0.0238 1818.3626 0.9993 ⋯
w[16] -2.7241 0.9117 0.0288 0.0175 2450.0885 0.9992 ⋯
a 104.2644 0.3407 0.0108 0.0100 952.1175 1.0038 ⋯
σ 6.0974 0.1558 0.0049 0.0038 1775.3299 0.9990 ⋯
1 column omitted
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
w[1] -2.5281 -1.2683 -0.7015 -0.1222 1.1072
w[2] -2.5470 -1.3879 -0.8532 -0.2840 0.8402
w[3] -0.6997 0.2898 0.8241 1.3810 2.3930
w[4] -0.6385 0.2026 0.6853 1.1511 2.0438
w[5] -1.0201 0.0450 0.6056 1.1317 2.2009
w[6] -2.6198 -1.7510 -1.2338 -0.7059 0.1801
w[7] -1.3066 -0.3101 0.2989 0.8075 1.7683
w[8] 0.1302 1.2344 1.8055 2.3861 3.4671
w[9] -1.6704 -0.5649 -0.0236 0.4951 1.4676
w[10] 0.4834 1.3700 1.9061 2.4198 3.3756
w[11] -0.8996 0.1316 0.6925 1.2162 2.2533
w[12] -0.1711 0.7374 1.2467 1.7692 2.7470
w[13] -0.1589 0.7075 1.2470 1.8217 2.8017
w[14] -2.2540 -1.2961 -0.7754 -0.1916 0.7812
w[15] -4.5643 -3.4532 -2.9137 -2.2586 -1.3194
w[16] -4.3992 -3.2736 -2.7774 -2.1793 -0.8404
a 103.5640 104.0427 104.2656 104.4941 104.9266
σ 5.8152 5.9908 6.0953 6.1964 6.4151